Coarse-grained Simulations of Chemical Oscillation in a Lattice Brusselator System 
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Accelerated coarse-graining (CG) algorithms for simulating heterogeneous chemical reactions on 
surface systems have recently gained much attention. In the present paper, we consider such an issue 
by investigating the oscillation behavior of a two-dimension (2D) lattice-gas Brusselator model. We 
have adopted a coarse-grained Kinetic Monte Carlo (CG-KMC) procedure, where m x m microscopic 
lattice sites are grouped together to form a CG cell, upon which CG processes take place with well- 
defined CG rates. We find that, however, such a CG approach almost fails if the CG rates are 
obtained by a simple local mean field (s-LMF) approximation, due to the ignorance of correlation 
among adjcent cells resulted from the trimolecular reaction in this nonlinear system. By properly 
incorporating such boundary effects, we thus introduce the so-called b-LMF CG approach. Extensive 
numerical simulations demonstrate that the 6-LMF method can reproduce the oscillation behavior 
of the system quite well, given that the diffusion constant is not too small. In addition, we find that 
the deviation from the KMC results reaches a nearly zero minimum level at an intermediate cell 
size, which lies in between the effective diffusion length and the minimal size required to sustain a 
well-defined temporal oscillation. 

PACS numbers: 05.10.Ln, 82.40.Bj, 02.70.Uu 



I. INTRODUCTION 



When driven far from thermal equilibrium, heteroge- 
neous surface chemical reaction systems often show a 
variety of complex dissipative structures such as oscilla- 
tions, Turing patterns, spiral waves and turbulence [ll-Ho| . 
Traditionally, these phenomena are mainly observed at 
macroscopic scales ranging from a few to several hundred 
micrometers, but recent studies showed that they are 
also present in nanoscale systems. At the present time, 
two different theoretical approaches are used to describe 
these nonlinear behaviors of surface reactions. Mean 
field deterministic equations (MDFE), like the reaction- 
diffusion equation, can provide good qualitative descrip- 
tion of spatiotemporal dvnamics[lll Il2j. However, they 
are essentially phenomenological and neglect microscopic 
mechanisms such as lateral interactions between adsor- 
bate molecules. In addition, MFDE also ignores the 
molecular fluctuations which may play important roles 
in mesoscopic systems. Another approach, microscopic 
lattice models, takes explicitly into account the adsorp- 
tion, desorption, diffusion and reaction processes as ran- 
dom events, and one can use kinetic Mente Carlo (KMC) 
methods to yield detailed valuable information about 
the microscopic reaction properties [l3rll8| . This micro- 
scopic method can account for the molecular interactions 
and fluctuations directly, however, memory and speed 
of available computers limit the maximal spatial size of 
the system, which renders the direct KMC simulation 
of nanoscale spatiotemporal structures and populations 
a difficult task. Therefore, a promising way is to de- 
velop coarse-grained (CG) approaches, bridging the gap 
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between those two, aiming at significantly reducing the 
degree of freedom to accelerate the simulation on large 
length scale while properly preserving the microscopic 
fluctuation information and correct dynamics. 

Very Recently, two kinds of CG methods based on 
lattice-gas model have been proposed. One is the con- 
tinuum mesoscopic modeling devoloped by Mikhailov 
et.al. [191 - 124} which is derived from coarse-graining of 
the underlying microscopic master equation to get the 
functional Fokker-Planck equation and its correspond- 
ing stochastic partial differential equations (SPDE). The 
other is a discrete CG approach proposed by Vlachos and 
coworkers poT432l | , which is a kind of CG-KMC algorithm 
by grouping the microscopic lattice sites into coarse cells 
and the CG system evolves by a sequence of CG events 
associated with the microscopic processes. By numeri- 
cally solving the SPDE, Mikhailov has successfully in- 
vestigated the nucleation of single reactive adsorbate in 
one-dimension (ID) and 2D systems [l9[, the formation 
of stationary microstructures for single species with at- 
tractive lateral interactions in a 2D system[2(| and the 
formation of traveling nanoscale structures in a model of 
two different species on 2D reaction surface [2l[ . Vlachos 
et.al had in-depthly investigated the validity of the CG- 
KMC approach in ID conceptional systems with different 
potential form (2||[27]], the ID Ising system with spin ex- 
change [28 1 . p rototype model of ID diffusion through a 
membrane[29(, steady pattern formation of simple reac- 
tion model on 2D surface (30|. and so on. Furthermore, 
the authors had also discussed the possibility for the ex- 
tension of CG-KMC to complex lattices, multicomponent 
systems [3l[ and heterogeneous plasma membranes f32j . 
It has been shown that both CG methods enable dynamic 
simulations over large length and time scales and can ac- 
curately capture transient and equilibrium solutions as 
well as noise properties especially for long-ranged poten- 
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tials. II. MODEL AND METHOD 

A. The Brusselator Model 



In the present paper, our mainly purpose is to find an 
effective CG-KMC method to investigate the nonlinear 
oscillation behaviors of surface chemical system which 
has already been developed into a field of very active 
research HH ■ Here, we adopt a 2D lattice-gas Brus- 
selator model, which is a typical oscillatory system with 
a nonlinear autocatalytic trimolecular reaction. To pre- 
serve the microscopic information correctly and acceler- 
ate the simulation at the same time, a CG-KMC algo- 
rithm by grouping m x m microscopic sites into a CG 
cell is adopted. Numerical results show that such a CG 
method is actually not a good approximation if the CG 
rates are obtained by a simple local mean field (s-LMF) 
approximation, due to the correlation among adjacent 
cells resulted from the trimolecular reaction cannot be 
neglected. We thus proposed a 6-LMF approach, which 
has properly accounted for the boundary corrections to 
the CG reaction rates. Extensive numerical simulations 
show that the 6-LMF method can reproduce quite well 
the oscillation behaviors, obtained from the KMC simu- 
lations on the microscopic lattice. To quantitatively in- 
vestigate the accuracy of the CG methods, we introduce 
deviation coefficients ja for the oscillation amplitude and 
7t for the oscillation period between the CG-KMC and 
KMC, respectively. We find again that the 6-LMF is 
remarkably better than s-LMF, and there is an interme- 
diate CG cell size where the deviation reaches a nearly 
zero minimal level. We suggest that for the CG method 
to work, the cell size should lie in between the effective 
diffusion length and the minimal size required to sustain 
a well-defined temporal oscillation. 



We consider a modified Brusselator model on the 2D 
surface lattice as follows: 
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Herein, the reactions under consideration are assumed 
to run on a N x N square lattice. Sites are either va- 
cant (denoted by * ) or occupied by single U or V parti- 
cles(See Fig. fa). The subscripts 'g' and 'a' represent the 
species in gas phase and adsorbed on the surface, respec- 
tively. Process (1) denotes the adsorption of species U, 
(2) the desorption of U, (3) the conversion from adsorbed 
U to V, and (4) the desorption of V, respectively. Step 
(5) is the autocatalytic reaction wherein an adsorbed V 
molecule with two nearest neighbor U molecules converts 
to U. The parameters K a (a = 1, ...,5) represent dimen- 
sionless rate constants. Steps (6) and (7) denote the dif- 
fusion processes of U and V, respectively, where Kq and 
K-j are the corresponding diffusion constants. In Table I, 
these processes and their corresponding propensity func- 
tions are listed. Herein, Wi a (a = 1,...,7) represent the 
propensity function of the a-th process taking place at 
site i. The occupation function af (where <j> = U or V) 
denotes the state of a given surface site i: af = I if site 
i is occupied by species (j> and otherwise. Note that a 
vacant site is necessary for the adsorption of U or the dif- 
fusion processes. In the reaction process (5), both sites j 
and k must be nearest neighbors of site i. 



The paper is organized as follows. In Section [Hi we 
present the lattice Brusselator model and describe the 
methods in detail, including the KMC and CG-KMC pro- 
cedures. The numerical results are given in lllll where we 
mainly focus on the comparison between KMC and CG- 
KMC, by investigating the deviations in the oscillation 
amplitude and period as functions of the control param- 
eter. We end by conclusions in IIVI 



B. MF Description 

Assuming that the surface is well-mixed by the diffu- 
sion process, one may describe the system dynamics by 
the following MF deterministic equations, 

du . dv 

— - = (wi - w 2 - w 3 + w 5 ), -j- = (w 3 - w A - w 5 ), (8) 
at at 

where u and v are respectively the surface coverage of 
species U and V. iu a =i,...,5 denote the rate of reaction-a 
with 

W\ = K\{1 — u — v), w 2 = K 2 u, w 3 = K 3 u, 

W4 = K4V,W5 — K4VU 2 . (9) 
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TABLE I: Stochastic processes and corresponding reaction rates in KMC simulation for the lattice-gas Brusselator model. 
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In the present work, we use this MF equation to deter- 
mine the bifurcation diagram of the system. In certain 
parameter ranges, the system can undergo a Hopf bifur- 
cation, where stable limit cycle emerges. 

The MF equation does not take into account inter- 
nal fluctuations inherent in chemical reaction systems. 
For small systems, such fluctuations may play important 
roles. To account for the internal noise while keep the 
MF approximation, one may use the chemical Langevin 
equations as follows, 

du 

— ={wi - w 2 -w 3 + w 5 ) 

+ jjWwiti® - (t) - v^3&(*) + VSfc&W 

(10) 

dv 

— ={w 3 -w A - w 5 ) 

+ jfWwfa{t)- V^Ut) - V^Ut)] (11) 

where £, a =i,...,5(t) are Gaussian white noises associated 
with the reaction channels, obeying < £ a (t) >= and 
< £a(t)£p(t') >= 5 a p5(t — t 1 ). The items in the bracket 
with £, a (t) give the internal noises, scaling as 1/N, which 
are closely coupled with the reaction kinetics. In the 
macroscopic limit N — > 00, the internal noise items can 
be ignored and the CLE recovers the deterministic equa- 
tion (8). 

C. KMC Simulations 

Given the processes and their propensity functions as 
listed in Table I, one can then perform KMC to study 
the dynamics. In the present work, we adopt a null- 
event KMC procedure. The main steps can be outlined 
as follows, 

1. Determine which process a to happen. To do this, 
we first draw a random number r\ from the uniform 
distribution in the unit interval, and then take a as 



the smallest integer satisfying ^3=1-^8 > r iW / o J 

where Wq = Y%—x Kp denotes the total maximum 
transition rates. 

2. Randomly select a surface site i with equal proba- 
bility pi = 1/N 2 . 

3. Determine whether the selected process a can take 
place on site i or not. This is given by a so-called 
participation index 34] €i a — W{ a j K a , which takes 
value or 1, corresponding to the propensity func- 
tions Wi a shown in Table I. Note that this index de- 
pends on the local configuration around site i and 
varies with time. For the desorption process (2), 
for instance, = 1 if site i is occupied by U and 
otherwise. For the diffusion of [/(or V), Ci a = 1 if 
the site i is occupied by U (or V) and a randomly 
select nearest neighbor j is vacant. For the reaction 
process (5), the index reads 1 if site i is V, and two 
succeedingly selected nearest neighbors j and k are 
both U. Note here we have used the same rule for 
this process as that proposed by Zhadnov [35|: As 
shown in Fig.l, reaction can only happen among 
orthogonal configurations, such as those shown by 
(a), (b) and (c), but not within line configurations 
such as (d). 

4. Execute the process a if e; Q = 1, and terminate the 
trial otherwise. 

5. Repeat steps 1 to 4. 

In the present work, we start the KMC run from a 
clean surface. The time advancement can be measured in 
terms of tkmc — 1/Wn. The KMC results are assumed 
to be correct, and we use them to check the validity of 
other methods, especially that of the CG-KMC. 

D. The CG-KMC Methods 

In the present work, we adopt a CG procedure origi- 
nally introduced by Vlachos et.al (28|. Since the diffusion 
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FIG. 1: The scheme of spatial coarse graining is plotted. Here, 
a 18 x 18 microscopic lattice (see dash-line) is divided into a 
3x3 CG lattice(see solid-line) by uniformly grouping 6x6 
microscopic sites to a CG cell. At the microscopic level, the 
reaction 2U + V — > 3U can only take place on orthogornal 
UVU configurations (a), (b) and (c), but not on line configu- 
ration (d). 



processes are usually faster than other slow processes, it 
is reasonable to assume that particles are well-mixed in 
a comparatively large domain, whose scale is determined 
by the diffusion length, over short time scales. There- 
fore, one can divide the micro-lattice into several coarse 
cells, wherein each particle is assumed to have an equal 
probability of occupying any microscopic lattice site, and 
no spatial correlation exists. Obviously, the most natural 
way for such a spatial coarse graining on 2D surface is 
to group to x to sites into a CG cell. For instance, Fig.l 
shows the coarse-graining of a 18 x 18 micro-lattice into a 
3x3 CG-lattice, wherein each CG-cell contains to 2 = 36 
sites (hereafter, we use 'site' for the micro-lattice and 
'cell' for the CG-lattice). 

To perform CG-KMC, one needs to define the CG pro- 
cesses taking place on the CG lattice and obtain the cor- 
responding CG rates. We now introduce the CG variables 

4=H 4 and 4 = 4/ m2 

to denote the number and coverage of ^-species in the [i- 
th CG-cell C M , respectively. Consider the micro- process 
(1), for example, the CG process is defined as the ad- 
sorption of one U particle in any site i inside a CG-cell 
Cfj,. Since [/-adsorption only involves one surface site, 
there is no spatial correlation between different adsorp- 
tion events inside a cell, therefore the rate of the CG 
adsorption process easily reads 



Following this simple rule, one can readily obtain the CG 
rates for the single-site processes (1) to (4). 

For the diffusion processes (6) and (7), one must bear 
in mind that the diffusion constant between two adja- 
cent CG-cell, Kq ( Kj ), is not identical to that be- 
tween two adjacent micro-sites, Kq{K-;). To establish 
the relationship between these two rates, we can adopt 
the so-called 'flux-consistency' rule, which requires that 
the average flux across the boundary of two adjacent CG 
cells, calculated from the CG diffusion process, should 
be the same as that calculated from averaging over the 
micro-diffusion process. By using this criterion and con- 
sidering the maintenance of detail-balance, see [ItJ, one 
must have Kq — K^/m 2 and K-j = K^/m 2 . Finally, the 
rate for the CG diffusion process of U is 

w, 6 =J2K 6 {aY(l-af-ay)) ieCM - ^<(1-^-^). 

i 

Herein, (•) means ensemble average. In the final equation, 
we have ignored the correlation between different cells 
and simply replaced the ensemble average of af inside 
C M by fj* |1 [i^. The rate for CG diffusion of V can be 
obtained in a similar way. 

For the trimolecular reaction process (5), however, 
strong correlation exists between neighboring sites. To 
perform the CG simulation, we need to express the 
ensemble-averaged rate of this process inside a cell by 
a function of the CG-variables, 

w m5 = F(K 5 ,m,f}V,fjV). 

However, it is hard to decide the correct functional form 
at this stage. It is worthy to note here that a seamless 
approach has been proposed very recently to address the 
validity of reaction-diffusion master equations 36] . The 
authors argued that to make the master equation to be 
consistent with the micro-model, the reaction constant 
must be dependent on the coarse-size, here is m. They 
demonstrated the success of this idea for a reversible 
aggregation-dissociation reaction. Unfortunately, it is 
rather difficult for us to work out a similar result for the 
trimolecular lattice gas system considered here. There- 
fore, to step forward, we have to make some approxima- 
tions. To the lowest order, one may use the simple local 
mean field (s-LMF) approximation as follows, 

UV5 = m 2 {K 5 ayaVaUj,kec, = m 2 K 5 ^(fj^) 2 . 

Herein, (•)' denotes the ensemble average over the orthog- 
onal configurations inside C M . Given that the diffusion 
process is fast and the cell size is smaller than the dif- 
fusion length, this approximation, although crude, might 
be acceptable. 

Unfortunately, as we will show(see below), however, 
this s-LMF scheme for the reaction process almost fails 
to reproduce the KMC results, no matter how large the 
coarse-size m and the diffusion constants Kq(K^) are. It 
seems that the s-LMF loses some key components that 
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should be considered during the CG procedure. We note 
here that the s-LMF scheme totally ignores the correla- 
tions between adjacent cells resulted from the trimolecu- 
lar process. For instance, the reaction configuration (b) 
(see in Fig.l) on the border of cell \i involves two sites 
(U and V) in cell (i and one site (U) in cell v. If we use 
LMF for both cell \x and v, the ensemble averaged rate 
for this particular configuration should read Ksffffj^fj^, 
which is different from the rate for the reaction inside C M , 
K 5 (fjU) 2 f]V ^ if W e consider that concentration gradients 
exist between adjacent cells. We argue that this effect 
should be taken into account to compensate the discrep- 
ancy resulted from the CG approximation. Similarly, the 
reaction configuration (c)(in Fig. 1) at the corner of cell 
H involves sites in three adjacent cells. Therefore, instead 
of the s-LMF, one may use a boundary-corrected LMF 
(b-LMF) scheme by writing down the CG-rate of process 
(5) as follows, 

corner- 

(12) 

Herein, 

w mtra ^m 2 K 5 fjl(fj^) 2 , (13) 

Wborder = m 2 K 6 fj^ (£, *^/ 4 )' ( 14 ) 

V 

and 

w corner = m 2 K 5 f)KJ2 (15) 

denote the average rate of process (5) inside, on the bor- 
der of, and at the corner of cell fi, respectively. Note that 
the summation in equation (14) runs over the four adja- 
cent cells of C M , and that in (15) runs over adjacent cells 
of the four corners. The weighting factors /i, / 2 and / 3 
denote the possibility of finding an reaction UVU con- 
figuration belonging to the three categories, respectively, 
given that the V site is inside the current cell C^. By 
simple manipulations, we have 

fi = (1 - I/™) 2 , h = 1/m 2 , and f 2 = 1 - h - f 3 . 

In Table II, the CG processes as well as their corre- 
sponding CG rates are listed. Note that b-LMF and s- 
LMF show difference only for the trimolecular reaction. 
According to these processes and rates, one can readily 
perform CG-KMC simulations. In the present paper, we 
also use null-event procedure as that for the KMC. The 
steps are outlined as follows, 

1. Choose a process similar to the first step used in 
the KMC, except that now Kq and K-j should be 
replaced by K 6 = K§/m 2 and K-j — K 7 /m 2 . Cor- 
respondingly, Wo should be changed to W GG . 

2. Randomly select a cell \i with equal probability. 



3. Calculate the reaction probability e MQ for the pro- 
cess a to happen associated with the current cell 
[i. This probability simply equals to w^/ (K a m 2 ) 
for 1 < a < 5 and w^j (K a m 2 ) for a = 6 or 7. 

4. Generate a second uniformly distributed random 
number r 2 in the unit interval. If r 2 < e MQ , execute 
the process a, and the trial ends otherwise. 

5. Repeat the above steps. 

In the present study, we start CG-KMC simulations 
from the same initial conditions as in the KMC. To be 
consistent with the KMC, the time increment should read 
t cg = 1/Wq G for each trial. We compare the CG-KMC 
results with the KMC ones to check the validity of CG 
approaches. 



III. NUMERICAL SIMULATIONS AND 
DISCUSSION 

In our work, the main parameters used in simulations 
are K x = 5.0 x 1(T 5 , K 2 = 1.0 x 10~ 3 , K 3 = 5.0 x 10~ 3 
and K 4 = 6.0 x 10~ 5 , while K 5 and K 6 = K 7 = D are 
control parameters. To compare the results of different 
methods, we have generated time series u(t) or v(t) with 
enough length and calculated the oscillation amplitude 
A and period T as a function of the control parameters. 
In Fig. (2a), the dependence of the oscillation range of v 
on K 5 is shown, obtained by the MFDE, CLE and KMC 
with different diffusion constant D. Correspondingly the 
curves for the period T are drawn in Fig. (2b). The CLE 
results are obtained by numerical simulation of Eq.(10) 
and (11) with a time step At = 0.01 and N = 256. 
All the KMC results are also performed on a 256 x 256 
square lattice. The solid line obtained from MFDE cor- 
responds to the bifurcation diagram. A Hopf bifurcation 
locates at K 5c ~ 2.15, below which deterministic oscil- 
lation can be observed. Several points can be addressed 
from this figure. First of all, CLE and KMC show strong 
qualitative differences with the MFDE: Stochastic oscil- 
lations can be observed even outside the deterministic 
oscillatory region, here K 5 > K 5c . This so-called noise 
induced oscillation phenomenon has gained great atten- 
tion in recent years and may have important applications 
especially in circadian oscillation systems. Secondly, the 
KMC results depend strongly on the diffusion constant 
D. However, the results for D = 10 and D = 30 nearly 
collapse, indicating that the KMC results may converge 
in the limit of large D. In this latter case, the MFDE can 
reproduce the KMC results when the parameter lies deep 
inside the oscillatory region, see the range K§ < 2.1. If D 
is small, both MFDE and CLE show large discrepancies 
with the KMC results, no matter the range of the control 
parameter. Finally, we would like to point out here that 
the CLE cannot reproduce the KMC results accurately 
even for large D, although they share some qualitative 
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TABLE II: CG processes and rates associated with the CG cell C M . 



Process Description State Change s-LMF Rate b-LMF Rate 

U Adsorption if -s> + 1 = m 2 K 1 (l - fj" — f)*) w^i 

U Desorption rff — ► 77]^ — 1 u> M 2 = rn^K^f^ mj m 2 

U Conversion to V ->■ 77^ — 1, 77^ -> jjjf + 1 «i m3 = m 2 K 3 fj^ uj m3 

V Desorption — s- ry^f — 1 uj m 4 = m 2 K4fjY «) M 4 

2U+V Reaction ^ -»• - 1, 77^ 77^ + 1 ™ m5 = m 2 K 5 fj^ (77^ ) 2 tu£ 6 

U Diffusion 77^ ->• r?^ - 1, 77^ -+rfi + 1 mj m6 = K e fj^ (1 - fff - ftf) uj m6 

V Diffusion v l ^1% - ljtf -> ^ + 1 Wp7 = jfrfff (1 - 77? - 77^) «i M 7 
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FIG. 2: The oscillation range of (a) V-coverage and (b) period 
are presented as functions of control parameter K^, obtained 
by MFDE, CLE and KMC with D = 1, 10, 30, respectively. 
Parameters are K x = 5.0 x 10~ 5 , K 2 = 1.0 x 10~ 3 , K 3 = 
5.0 x 10~ 3 , K 4 = 6.0 x 10~ 5 and N = 256. 



FIG. 4: The oscillation range of (a) V-coverage and (b) pe- 
riod obtained from the 6-LMF CG approach with different 
coarse size m. The diffusion constant is D = 10 and the other 
parameters are the same as in Fig. 2. 




FIG. 3: The oscillation range of (a) V-coverage and (b) pe- 
riod obtained from the s-LMF CG approach with different 
coarse size m. The diffusion constant is D = 10 and the other 
parameters are the same as in Fig. 2. 



features, e.g., noise induced oscillation to the right side 
of the Hopf point. 

In the following part, we mainly consider a system 
with size N = 256 and diffusion constant D = 10. We 
have also performed some KMC simulations on larger 
systems, e.g., N = 512, but the main conclusions are the 



same. Since we are mainly interested in the validity of 
CG methods and extensive simulations are required to 
compare the results of different methods, we have fixed 
N = 256 throughout the paper. CG-KMC simulations 
are performed according to the CG processes and rates 
listed in Table II. In Fig. 3, the oscillation amplitude and 
period obtained by using s-LMF rates are shown, for dif- 
ferent sizes of the CG cell. Apparently, the s-LMF ap- 
proach almost fails to reproduce the results of KMC, even 
qualitatively. For small to, the s-LMF totally loses the 
whole bifurcation features of the KMC dynamics. We 
show in Fig. 4, however, that the b-LMF behaves much 
better than the s-LMF. Firstly, the b-LMF can repro- 
duce the global bifurcation feature quite well, even for 
small to. In addition, for an intermediate value of to, say, 
to = 8 here, the b-LMF results show excellent consistent 
with the KMC results, in both the oscillation amplitude 
(Fig. 4a) and the period (Fig. 4b). It seems that the b- 
LMF approach does catch some key factors during the 
CG procedure. 

In Fig. 5, we have plotted the dependence of the 
turnover frequency (TOF) as a function of time for dif- 
ferent coarse size to obtained from the KMC, b-LMF 
and s-LMF. The TOF is denned as the occurrence of the 
trimolecular reaction (5) per surface site per unit time. 
Clearly, the b-LMF with to = 8 matches the KMC quite 
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FIG. 6: The dependence of the deviation coefficient (a) for 
the oscillation amplitude 7a and (b) oscillation period 7t on 
the coarse cell size m obtained from the s-LMF and 6-LMF 
methods. Parameter are the same as in Fig. 2 and N = 256. 



well, while that with smaller or larger m may capture 
some qualitative features of the TOF but with apparent 
quantitative differences. The s-LMF, however, almost 
loses the temporal information associated with the TOF. 

To further demonstrate this quantitatively, we intro- 
duce a deviation coefficient for the oscillation range as 
follows. According to Fig. 3 and Fig. 4, each bifurcation 
diagram contains two branches, the upper branch and 
the lower one corresponding to the averaged maximum 
and minimum values of v(t), respectively. As can be seen 
from the figures, both branches obtained from the CG- 
KMC methods show discrepancies with the KMC values. 
Denote the upper branch value of v, obtained by CG- 
KMC, at a certain control parameter if 5 by v%, and that 
obtained by KMC by v% , then 



measures the relative discrepancy of the upper branch, 
where Nk is the number of control parameters used in the 
calculation. Similarly, we can calculate the discrepancy 
of the lower branch j l A . In the present work, we have 
used Nk = 20 points inside the range K 5 e (1.75,2.5) 
to obtain 7^ and j l A . In Fig. (6a), the dependence of 
ja = (Ja + Ta)/^ on the size m of CG-cell is shown, for 
the b-LMF with different diffusion constant D and the 
s-LMF with D = 10. In Fig. (6b), the curves for 77-, the 
relative discrepancy in the period, are shown. Clearly, 
the s-LMF method shows relatively large discrepancies, 
while the b-LMF works much better. The s-LMF is even 
worse than the MFDE, shown by the dash lines in Fig. 6 
for D = 10. One notes that both ja and jt exhibit a 
clear-cut minimum of about zero at m = 8 for large D 
when b-LMF is used. We also note that if D is small, say 
D = 1 here, the b-LMF also fails. This is not surprising 
because CG method which assumes well-mixing in a CG- 
cell should not work if diffusion is too slow. 



In the above results, we see that the CG results for 
small m do not match the results of KMC, even for the 
b-LMF. This is in contrast to the CG-KMC methods used 
by Vlachos et.al to account for the dynamics of 2D lat- 
tice gas Ising model. Note that for the Ising model, one 
mainly considered the equilibrium states. For the Brus- 
selator model considered here, however, we want to re- 
produce the temporal oscillation behavior. To reproduce 
the oscillation features on the whole surface, the tempo- 
ral correlation of the time series must be properly main- 
tained during the CG procedure. When we perform the 
CG procedure by dividing the lattice into CG cells, we 
are dealing with A^ c x N c coupled CG oscillators, where 
N c = N/m is the size of the CG lattice. The coupling 
between these CG oscillators are realized by the diffusion 
and the boundary correlation considered in the b-LMF. 
If the CG cell is too small, however, the time-correlation 
inside each cell will be lost due to strong fluctuations and 
the time evolution of 77^ (t) cannot be viewed as an oscil- 
lation. As discussed by P.Gaspard, a minimum number 
of well-mixed molecules is required to produce correlated 
oscillations, such that the auto-correlation time of the 
time series is not smaller than T/27r[3^, Therefore, 
it seems that a seamless CG approach to reproduce tem- 
poral dissipative structures like chemical oscillation is a 
large challenge. On the other hand, for any CG method 
within LMF scheme to work well, the scale of a CG cell 
should not be larger than the diffusion length, as empha- 
sized by Mikhailov and others 0, The compromise 
between these two factors, i.e., to keep time autocorrela- 
tion and to be smaller than the diffusion length, may be 
the reason of appearance of an optimal m for the b-LMF 
approach. We note that this reasoning is not applicable 
to the s-LMF, for which the discrepancies monotonically 
decrease with increasing m, since the s-LMF does not 
work for the present system. 
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IV. CONCLUSIONS 

In summary, we have tried to apply a CG-KMC ap- 
proach to simulate the oscillation behavior on the sur- 
face lattice-gas Brusselator model. Owing to the correla- 
tions between adjacent cells resulted from the nonlinear 
trimolecular reaction, the CG approach based on sim- 
ple LMF approximation almost fails. By properly taking 
into account the boundary corrections, we have intro- 
duced a so-called b-LMF CG-KMC approach, which can 
reproduce the microscopic KMC results quite well, given 
that the diffusion is not too slow and the CG cell size is 



optimally chosen. Our work thus unravels the very role 
of reaction correlations which should be carefully consid- 
ered in any CG approach and mesoscopic modeling for 
nonequilibrium spatiotemporal dynamics at nanoscales. 
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